Prostorni podaci u R-u - vežbe (drugi deo)

Uvod

U okviru ovih vežbi biće predstavljene mogućnosti definisanja i manipulacije nad definicijom referentnog koordinatnog sistema prostornih podataka u R-u. Zatim biće predstavljene mogućnosti učitavanja, obrade, analize i manipulacije nad rasterskim podacima.

Koordinatni referentni sistem

Coordinate reference system (CRS) - Koordinatni referentni sistem ili Prostorni referentni sistem - Spatial reference system (SRS) predestavlja sinonim kojim se definiše niz parametara kojim se definiše položaj određene tačke u prostoru u odnosu na položaj Zemlje.

Geografski koordinatni referentni sistem je geometrijski model u kome je definisano sledeće:

1. Model oblika zemlje (npr: elipsoid sa parametrima \(a\) i \(e\));
2. Početni meridijan (najčešće Grinički meridijan, pa se ovaj deo podrazumeva i često izostavlja iz definicije);
3. Datum.

Koordinatni referentni sistem u projekciji je geometrijski model u kome je definisano sledeće:

1. Projekcija;
2. Model oblika zemlje (npr: elipsoid sa parametrima \(a\) i \(e\));
3. Početni meridijan (npr: Grinički meridijan);
4. Datum.

Prema Lott-u 2015 važno je razlikovati:

• Koordinatni sistem (coordinate system) je set matematickih pravila koja definišu kako se koordinate dodeljuju tacki
• Datum (datum) je set parametara koji definiše položaj centra, razmeru i orijentaciju koordinatnog sistema
• Geodetski datum (Geodetic datum) opisuje odnos koordinatnog sistema u odnosu na Zemlju
• Koordinatni Referentni Sistem (Coordinate Reference System) je koordinatni sistem ciji je položaj odreden datumom u odnosu na referentni objekat. U slucaju planete zemlje (WGS84) parametri geodetskog datuma se koriste.
• Možemo razlikovati Geografske koordinatne sisteme, Sisteme u projekciji i Lokalne koordinatne sisteme.

Učitavanje potrebnih paketa: ***

library(sp)
library(raster)
library(sf)
library(mapview)
library(stringr)
library(tidyverse)
library(dplyr)
library(magrittr)
library(RStoolbox)
library(gridExtra)
library(viridis)
library(stars)

Učitavanje podataka

---

gradovi <- st_read(system.file("vectors/cities.shp", package = "rgdal"))

## Reading layer `cities' from data source `C:\Users\pbursac\Documents\R\win-library\4.0\rgdal\vectors\cities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 606 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84

gradovi

## Simple feature collection with 606 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84
## First 10 features:
##                NAME COUNTRY POPULATION CAPITAL                   geometry
## 1          Murmansk  Russia     468000       N  POINT (33.08604 68.96355)
## 2       Arkhangelsk  Russia     416000       N  POINT (40.64616 64.52067)
## 3  Saint Petersburg  Russia    5825000       N  POINT (30.45333 59.95189)
## 4           Magadan  Russia     152000       N      POINT (150.78 59.571)
## 5             Perm'  Russia    1160000       N  POINT (56.23246 58.00024)
## 6     Yekaterinburg  Russia    1620000       N  POINT (60.61013 56.84654)
## 7  Nizhniy Novgorod  Russia    2025000       N  POINT (43.94067 56.28968)
## 8           Glasgow      UK    1800000       N POINT (-4.269948 55.86281)
## 9            Kazan'  Russia    1140000       N  POINT (49.14547 55.73301)
## 10      Chelyabinsk  Russia    1325000       N    POINT (61.39261 55.145)

Vizuelizacija podataka

---

mapview(gradovi, zcol = "POPULATION")

Kreiaranje data frame-a gde su prostorni podaci - koordinate tačaka date kao dve kolone:

gradovi_df <- cbind(st_drop_geometry(gradovi), 
                    st_coordinates(gradovi))

head(gradovi_df, 10)

##                NAME COUNTRY POPULATION CAPITAL          X        Y
## 1          Murmansk  Russia     468000       N  33.086040 68.96355
## 2       Arkhangelsk  Russia     416000       N  40.646160 64.52067
## 3  Saint Petersburg  Russia    5825000       N  30.453327 59.95189
## 4           Magadan  Russia     152000       N 150.780014 59.57100
## 5             Perm'  Russia    1160000       N  56.232464 58.00024
## 6     Yekaterinburg  Russia    1620000       N  60.610130 56.84654
## 7  Nizhniy Novgorod  Russia    2025000       N  43.940670 56.28968
## 8           Glasgow      UK    1800000       N  -4.269948 55.86281
## 9            Kazan'  Russia    1140000       N  49.145466 55.73301
## 10      Chelyabinsk  Russia    1325000       N  61.392612 55.14500

sf paket - definicija CRS-a

---

Definisanje CRS-a se postiže na sledeći način:

crs_wgs84 <- st_crs(4326) # WGS84 elipsoidni koordinatni sistem ima EPSG kod 4326
class(crs_wgs84) # klasa

## [1] "crs"

crs_wgs84 # kod crs objekta uočavaju se sve stvari: EPSG kod i WKT2 string

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

cat(crs_wgs84$wkt) # direktan pristup wkt elementu crs objekta

## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

crs_wgs84$epsg # direktan pristup epsg elementu crs objekta

## [1] 4326

crs_wgs84$proj4string # direktan pristup proj4string elementu crs objekta

## [1] "+proj=longlat +datum=WGS84 +no_defs"

Dodeljivanje CRS-a se postiže na sledeći način:

gradovi2 <- st_as_sf(gradovi_df, coords = c("X", "Y")) # kreiranje objekta sf klase
gradovi2 # CRS: NA

## Simple feature collection with 606 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## CRS:            NA
## First 10 features:
##                NAME COUNTRY POPULATION CAPITAL                   geometry
## 1          Murmansk  Russia     468000       N  POINT (33.08604 68.96355)
## 2       Arkhangelsk  Russia     416000       N  POINT (40.64616 64.52067)
## 3  Saint Petersburg  Russia    5825000       N  POINT (30.45333 59.95189)
## 4           Magadan  Russia     152000       N      POINT (150.78 59.571)
## 5             Perm'  Russia    1160000       N  POINT (56.23246 58.00024)
## 6     Yekaterinburg  Russia    1620000       N  POINT (60.61013 56.84654)
## 7  Nizhniy Novgorod  Russia    2025000       N  POINT (43.94067 56.28968)
## 8           Glasgow      UK    1800000       N POINT (-4.269948 55.86281)
## 9            Kazan'  Russia    1140000       N  POINT (49.14547 55.73301)
## 10      Chelyabinsk  Russia    1325000       N    POINT (61.39261 55.145)

st_crs(gradovi2) <- 4326 # Putem EPSG koda!
# Drugi način: st_crs(gradovi2) <- crs_wgs84

# Pristup CRS objektu
st_crs(gradovi2)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

# Pristup wkt elementu CRS objekta
cat(st_crs(gradovi2)$wkt)

## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

sp paket - definicija CRS-a

---

Definisanje CRS-a se postiže na sledeći način:

crs_wgs84 <- CRS(SRS_string = "EPSG:4326") # WGS84 elipsoidni koordinatni sistem ima EPSG kod 4326
class(crs_wgs84) # klasa

## [1] "CRS"
## attr(,"package")
## [1] "sp"

crs_wgs84 # defnisan putem proj4string-a

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

wkt_wgs84 <- wkt(crs_wgs84)
cat(wkt_wgs84)

## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

Dodeljivanje CRS-a se postiže na sledeći način:

gradovi2 <- gradovi_df # kreiranje objekta sp klase
coordinates(gradovi2) <-  ~ X + Y
proj4string(gradovi2) <- crs_wgs84

gradovi2

## class       : SpatialPointsDataFrame 
## features    : 606 
## extent      : -165.27, 177.1302, -53.15, 78.2  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 4
## names       :    NAME,     COUNTRY, POPULATION, CAPITAL 
## min values  : Abidjan, Afghanistan,        -99,       N 
## max values  :    Zibo,    Zimbabwe,   23620000,       Y

str(gradovi2)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  606 obs. of  4 variables:
##   .. ..$ NAME      : chr [1:606] "Murmansk" "Arkhangelsk" "Saint Petersburg" "Magadan" ...
##   .. ..$ COUNTRY   : chr [1:606] "Russia" "Russia" "Russia" "Russia" ...
##   .. ..$ POPULATION: num [1:606] 468000 416000 5825000 152000 1160000 ...
##   .. ..$ CAPITAL   : chr [1:606] "N" "N" "N" "N" ...
##   ..@ coords.nrs : int [1:2] 5 6
##   ..@ coords     : num [1:606, 1:2] 33.1 40.6 30.5 150.8 56.2 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:606] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "X" "Y"
##   ..@ bbox       : num [1:2, 1:2] -165.3 -53.2 177.1 78.2
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "X" "Y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"

# Pristup CRS objektu u formi wkt elementa
cat(wkt(gradovi2))

## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

Konverzija sf <–> sp

---

gradovi_sf <- st_as_sf(gradovi2) # iz sp klase u sf klasu
gradovi_sf

## Simple feature collection with 606 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## geographic CRS: WGS 84 (with axis order normalized for visualization)
## First 10 features:
##                NAME COUNTRY POPULATION CAPITAL                   geometry
## 1          Murmansk  Russia     468000       N  POINT (33.08604 68.96355)
## 2       Arkhangelsk  Russia     416000       N  POINT (40.64616 64.52067)
## 3  Saint Petersburg  Russia    5825000       N  POINT (30.45333 59.95189)
## 4           Magadan  Russia     152000       N      POINT (150.78 59.571)
## 5             Perm'  Russia    1160000       N  POINT (56.23246 58.00024)
## 6     Yekaterinburg  Russia    1620000       N  POINT (60.61013 56.84654)
## 7  Nizhniy Novgorod  Russia    2025000       N  POINT (43.94067 56.28968)
## 8           Glasgow      UK    1800000       N POINT (-4.269948 55.86281)
## 9            Kazan'  Russia    1140000       N  POINT (49.14547 55.73301)
## 10      Chelyabinsk  Russia    1325000       N    POINT (61.39261 55.145)

gradovi_sp <- as(gradovi_sf, "Spatial") # iz sf klase u sp klasu
gradovi_sp

## class       : SpatialPointsDataFrame 
## features    : 606 
## extent      : -165.27, 177.1302, -53.15, 78.2  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 4
## names       :    NAME,     COUNTRY, POPULATION, CAPITAL 
## min values  : Abidjan, Afghanistan,        -99,       N 
## max values  :    Zibo,    Zimbabwe,   23620000,       Y

Prostorni podaci u formi raster-a

---

Raster ili grid je prostorna struktura podataka koja deli prostor u pravougaone ćelije (piksele) koje sadrže jednu ili više vrednosti nekog opažanog ili modeliranog fenomena, vrednosti ćelija su uglavnom numerički ili kategorijski podaci.
Više na: Geovizualizacija i Web kartografija, M. Kilibarda i D. Protić - http://osgl.grf.bg.ac.rs/books/gvvk/.

Raster (grid) koji ima jednu atributnu vrednost podseća na matricu koja dodatno ima definisano i zaglavlje. Postavlja se pitanje kako se definiše georeferencija kod rasterskih podataka. Praktično ona je sadržana u zaglavlju. Zavisno od formata rasterskih podataka, u zaglavlju se najčešće nalaze podaci o koordinatnom referentnom sistemu (kod nekih formata su u zasebnom fajlu), broju vrsta i kolonan matrice podataka, koordinatama početne ćelije (piksela), veličini piksela (neki formati podržavaju samo kvadratne dok drugi podržavaju i pravougaone piksele), informacija kako je zapisana vrednost za piksele koji nemaju atributnu vrednost.

U okviru R okruženja, za rad sa prostornim podacima u formi rastera, dostupni su paketi raster i stars.

Paket raster

Paket raster sadrži više klasa za rasterske podatke od kojih su tri glavne RasterLayer, RasterBrick i RasterStack. Raster paket sadrži funkcionalnosti za kreiranje, čitanje, pisanje i manipulaciju nad rasterskim podacima, kao što je na primer rasterska algebra. Isto tako paket sadrži i funkcionalnosti za vizualizaciju rasterskih podataka.

Klase:

• Klasa RasterLayer predstavlja rasterski podatak sa jednim atributom. Objekat tipa RasterLayer sadrži pored atributnih podataka i sve parametre koje se tipično čuvaju u zaglavlju (hederu) rasterskog fajla. To su: broj vrsta i kolona, prostorni obuhvat podataka i parametri koordinatnog referentnog sistema. Zavisno od veličine fajla svi podaci mogu biti u RAM memoriji ili im se pristupa sa diska računara.
• RasterBrick je klasa koja podržava više atributa u okviru rasterskih podataka, tipično to važi kad se koristi satelitski snimak sa više kanala, na primer multispektralni snimak iz Sentinel 2 misije ili Landsat 8. Ono što je specifično za RasterBrick je da on može da se referiše samo na jedan fajl koji je sačuvan na disku, na primer GeoTIFF fajl u kome su sačuvani kanali nekog multispektralnog satelitskog snimka.
• RasterStack je lista rasterskih lejera, vrlo slične strukture kao RasterBrick, ali ne postoji ograničinje da se referiše na jedan GIS fajl, već može da bude kompozicija više fajlova koji su sačuvani na disku ili koji su u RAM memoriji računara, ili kombinacija istih. Praktično se može reći da je to lista objekata tipa RasterLayer. Podrazumeva se da je prostorni obuhvat i rezolucija ista za svaki RasterLayer koji je deo RasterStack objekta.

Paket stars

Paket stars [https://r-spatial.github.io/stars/] nudi klase i metode za učitavanje, obradu, vizualizaciju i čuvanje podataka u formi vektorskih i rasterskih data cube - kocki podataka. U potpunosti se oslanja na tidyverse familiju paketa i logiku rada sa prostornim podacima koja je implementirana u okviru paketa sf.

Definisanje raster objekta i koordinatnog sistema

Kreiranje rastera sa slucajnim vrednostima celija - piksela

extent(7440500, 7475000, 4948000, 4976500) # xmin, xmax, ymin, ymax

## class      : Extent 
## xmin       : 7440500 
## xmax       : 7475000 
## ymin       : 4948000 
## ymax       : 4976500

podrucje_beograda <- raster(extent(7440500, 7475000, 4948000, 4976500), res = 100)

Funkcija rnorm - vektor slucajnih brojeva sa normalnom raspodelom

values(podrucje_beograda) <- rnorm(ncell(podrucje_beograda)) 
summary(values(podrucje_beograda))

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -4.037175 -0.675902 -0.000835 -0.000604  0.674081  4.119823

podrucje_beograda # CRS:NA

## class      : RasterLayer 
## dimensions : 285, 345, 98325  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 7440500, 7475000, 4948000, 4976500  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : -4.037175, 4.119823  (min, max)

plot(podrucje_beograda)

Definicija referentnog koordinatnog sistema

podrucje_beograda1 <- podrucje_beograda
podrucje_beograda2 <- podrucje_beograda
podrucje_beograda3 <- podrucje_beograda
podrucje_beograda4 <- podrucje_beograda

Dodeljivanje CRS-a se vrsi putem funkcije <- crs(), na vise načina koji daju isti rezultat

crs(podrucje_beograda1) <- 3909
crs(podrucje_beograda2) <- "EPSG:3909"
crs(podrucje_beograda3) <- st_crs(3909)$wkt # a WKT string
crs(podrucje_beograda4) <- CRS(SRS_string = "EPSG:3909") # an sp CRS object

Kontrola:

cat(wkt(podrucje_beograda1))

## PROJCRS["MGI 1901 / Balkans zone 7 (with axis order normalized for visualization)",
##     BASEGEOGCRS["MGI 1901",
##         DATUM["MGI 1901",
##             ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",3906]],
##     CONVERSION["Balkans zone 7",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",21,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9999,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",7500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",18277]],
##     CS[Cartesian,2],
##         AXIS["easting (Y)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["northing (X)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

all.equal(wkt(podrucje_beograda1),
          wkt(podrucje_beograda2),
          wkt(podrucje_beograda3),
          wkt(podrucje_beograda4))

## [1] TRUE

Vizualizacija:

podrucje_beograda1

## class      : RasterLayer 
## dimensions : 285, 345, 98325  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 7440500, 7475000, 4948000, 4976500  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -4.037175, 4.119823  (min, max)

mapview::mapview(podrucje_beograda1, na.color = "transparent", trim = TRUE) # trim - should rasters be trimmed off NA values around the edges

Multispektralni satelitski podaci

Dati su Multispektralni satelitski podaci - snimci, satelitske misije Sentinel 2, sa sledećim dostupnim kanalima:

U okviru ovog dela vežbe potrebno je:

• Učitati 3 multispektralna satelitska snimka (3 epohe), sa više dostupnih kanala.
• Kreirati raster stack kao kombinaciju sve tri epohe.
• Vizualizovati podatke u formi kolor kompozita i falš kolor kompozita.
• Sračunati vegetacioni indeks za sve tri epohe na nivou piksela.
• Kreirati histogram NDVI indeksa.

Učitavanje podataka

files <- list.files("Data/Multispektralni_snimci/", 
                    full.names = TRUE,
                    pattern = ".tif")

files

## [1] "Data/Multispektralni_snimci/0_2015-08-24_S2.tif"
## [2] "Data/Multispektralni_snimci/0_2015-09-13_S2.tif"
## [3] "Data/Multispektralni_snimci/0_2015-10-23_S2.tif"

Prvi način: svaki pojedinačan snimak kao poseban rasterStack putem for petlje

rasList <- list()
for(i in 1:length(files)){
  if(file.size(files[i]) > 0){
    rasList[[i]] <- raster::stack(files[i]) 
  } else {
    rasList[[i]] <- NULL
  }
}
rasList

## [[1]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.08.24_S2.1, X0_2015.08.24_S2.2, X0_2015.08.24_S2.3, X0_2015.08.24_S2.4, X0_2015.08.24_S2.5, X0_2015.08.24_S2.6, X0_2015.08.24_S2.7, X0_2015.08.24_S2.8, X0_2015.08.24_S2.9, X0_2015.08.24_S2.10, X0_2015.08.24_S2.11, X0_2015.08.24_S2.12, X0_2015.08.24_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535 
## 
## 
## [[2]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.09.13_S2.1, X0_2015.09.13_S2.2, X0_2015.09.13_S2.3, X0_2015.09.13_S2.4, X0_2015.09.13_S2.5, X0_2015.09.13_S2.6, X0_2015.09.13_S2.7, X0_2015.09.13_S2.8, X0_2015.09.13_S2.9, X0_2015.09.13_S2.10, X0_2015.09.13_S2.11, X0_2015.09.13_S2.12, X0_2015.09.13_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535 
## 
## 
## [[3]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.10.23_S2.1, X0_2015.10.23_S2.2, X0_2015.10.23_S2.3, X0_2015.10.23_S2.4, X0_2015.10.23_S2.5, X0_2015.10.23_S2.6, X0_2015.10.23_S2.7, X0_2015.10.23_S2.8, X0_2015.10.23_S2.9, X0_2015.10.23_S2.10, X0_2015.10.23_S2.11, X0_2015.10.23_S2.12, X0_2015.10.23_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535

Drugi način: svaki pojedinačan snimak kao poseban rasterStack putem lapply funkcije

rasStack <- lapply(files, stack)
rasStack

## [[1]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.08.24_S2.1, X0_2015.08.24_S2.2, X0_2015.08.24_S2.3, X0_2015.08.24_S2.4, X0_2015.08.24_S2.5, X0_2015.08.24_S2.6, X0_2015.08.24_S2.7, X0_2015.08.24_S2.8, X0_2015.08.24_S2.9, X0_2015.08.24_S2.10, X0_2015.08.24_S2.11, X0_2015.08.24_S2.12, X0_2015.08.24_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535 
## 
## 
## [[2]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.09.13_S2.1, X0_2015.09.13_S2.2, X0_2015.09.13_S2.3, X0_2015.09.13_S2.4, X0_2015.09.13_S2.5, X0_2015.09.13_S2.6, X0_2015.09.13_S2.7, X0_2015.09.13_S2.8, X0_2015.09.13_S2.9, X0_2015.09.13_S2.10, X0_2015.09.13_S2.11, X0_2015.09.13_S2.12, X0_2015.09.13_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535 
## 
## 
## [[3]]
## class      : RasterStack 
## dimensions : 28, 57, 1596, 13  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : X0_2015.10.23_S2.1, X0_2015.10.23_S2.2, X0_2015.10.23_S2.3, X0_2015.10.23_S2.4, X0_2015.10.23_S2.5, X0_2015.10.23_S2.6, X0_2015.10.23_S2.7, X0_2015.10.23_S2.8, X0_2015.10.23_S2.9, X0_2015.10.23_S2.10, X0_2015.10.23_S2.11, X0_2015.10.23_S2.12, X0_2015.10.23_S2.13 
## min values :                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                  0,                   0,                   0,                   0,                   0 
## max values :              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,              65535,               65535,               65535,               65535,               65535

Koordinatni sistem:

cat(wkt(rasStack[[1]]))

## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 34N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",21,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16034]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

for(i in 1:length(rasStack)){
  crs(rasStack[[i]]) <- 32634
}

cat(wkt(rasStack[[1]]))

## PROJCRS["WGS 84 / UTM zone 34N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 34N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",21,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World - N hemisphere - 18°E to 24°E - by country"],
##         BBOX[0,18,84,24]],
##     ID["EPSG",32634]]

Vizualizacija multispektralnih snimaka

---

Funkcija za vizualizaciju podataka u formi kolor kompozita i falš kolor kompozita:

Nazivi snimka - pomoćno:

files.names <- list.files("Data/Multispektralni_snimci/", full.names = F,
                          pattern = ".tif")
files.names %<>% stringr::str_remove(., pattern = ".tif")
files.names

## [1] "0_2015-08-24_S2" "0_2015-09-13_S2" "0_2015-10-23_S2"

plotS2parcel <- function(listS2parc = listS2parc, files.names = files.names, id = 1){
  
  c.plot <- ggRGB(listS2parc[[id]], r = 3, g = 2, b = 1) +
    labs(x = "Easting [m]", y = "Northing [m]",
         title = "Sentinel 2 satellite image - Color composite",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank())
  
  f.plot <- ggRGB(listS2parc[[id]], r = 7, g = 3, b = 2) +
    labs(x = "Easting [m]", y = "Northing [m]",
         title = "Sentinel 2 satellite image - Falsh color composite",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank())
  

  
  gridCF <- grid.arrange(c.plot, f.plot, ncol = 1, nrow = 2)
  return(gridCF)
}

Vizualizacija:

id1 <- plotS2parcel(listS2parc = rasStack, files.names = files.names, id = 1)

id2 <- plotS2parcel(listS2parc = rasStack, files.names = files.names, id = 2)

id3 <- plotS2parcel(listS2parc = rasStack, files.names = files.names, id = 3)

Računanje vegetacionog indeksa - NDVI

---

Pogledati: https://www.indexdatabase.de/

(NIR - Red) / (NIR + Red) Kanali: NIR = B08, i = 7 i RED = BO2, i = 3;

Funkcija za računanje NDVI indeksa nad svim snimcima:

rasListNDVI <- list()
for(i in 1:length(rasStack)){
    rasListNDVI[[i]] <- (rasStack[[i]][[7]] - rasStack[[i]][[3]]) / (rasStack[[i]][[7]] + rasStack[[i]][[3]])
}

rasListNDVI

## [[1]]
## class      : RasterLayer 
## dimensions : 28, 57, 1596  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0.1295294, 0.7217312  (min, max)
## 
## 
## [[2]]
## class      : RasterLayer 
## dimensions : 28, 57, 1596  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0.1917172, 0.7256516  (min, max)
## 
## 
## [[3]]
## class      : RasterLayer 
## dimensions : 28, 57, 1596  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 397599.5, 398169.5, 4987613, 4987893  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0.1663763, 0.8106579  (min, max)

summary(rasListNDVI[[1]])

##             layer
## Min.    0.1295294
## 1st Qu. 0.2300365
## Median  0.2941529
## 3rd Qu. 0.3932392
## Max.    0.7217312
## NA's    0.0000000

summary(rasListNDVI[[2]])

##             layer
## Min.    0.1917172
## 1st Qu. 0.2799979
## Median  0.4163690
## 3rd Qu. 0.5377995
## Max.    0.7256516
## NA's    0.0000000

summary(rasListNDVI[[3]])

##             layer
## Min.    0.1663763
## 1st Qu. 0.2133507
## Median  0.2371211
## 3rd Qu. 0.3487042
## Max.    0.8106579
## NA's    0.0000000

Histogram i vizuelizacija

---

plotS2parcel <- function(listS2parc = listS2parc, listS2parcNDVI = listS2parcNDVI, files.names = files.names, id = 1){
  c.plot <- ggRGB(listS2parc[[id]], r = 3, g = 2, b = 1) +
    labs(x = "Easting [m]", y = "Northing [m]",
         title = "Sentinel 2 satellite image - Color composite",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank())
    
  f.plot <- ggRGB(listS2parc[[id]], r = 7, g = 3, b = 2) +
    labs(x = "Easting [m]", y = "Northing [m]",
         title = "Sentinel 2 satellite image - Falsh color composite",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank())
  
  star_pred <- st_as_stars(listS2parcNDVI[[id]])
  
  ndvi.plot <- ggplot()+
    geom_stars(data = star_pred)+
    scale_fill_viridis(option = "D", na.value = NA, name = "NDVI: ")+
    labs(x = "Easting [m]", y = "Northing [m]",
         title = "Sentinel 2 satellite image - NDVI",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank(), 
          legend.position = "bottom")
  
  df <- as.data.frame(listS2parcNDVI[[id]])
  hist.plot <- ggplot(data=df, aes(layer)) + 
    geom_histogram(binwidth = 0.005,
                   col="red", 
                   aes(y = ..density..,
                       fill = ..count..)) +
    scale_fill_gradient("Count", low = "blue", high = "red")+
    stat_function(fun = dnorm, 
                  color = "black",
                  size = 1.5,
                  args = list(mean = mean(df$layer), sd = sd(df$layer)))+
    labs(x = "NDVI values", y = "Count",
         title = "Sentinel 2 satellite image - NDVI histogram",
         subtitle = paste("Image: ", files.names[id]),
         caption = "UBFCE (2020) \n
       Spatial resolution 10mx10m, Teritory of the Republic of Serbia") +
    theme(line = element_blank(),
          panel.background = element_blank(), 
          legend.position = "bottom")+
    annotate("text",x = 0.5, y = 5.5, label = paste("Mean:",round(mean(df$layer),digits = 2)))+
    annotate("text",x = 0.5, y = 5, label = paste("SD:",round(sd(df$layer),digits = 2)))
  
  gridCF <- grid.arrange(c.plot, f.plot, ndvi.plot, hist.plot, ncol = 2)
  return(gridCF)
}

id1 <- plotS2parcel(listS2parc = rasStack, listS2parcNDVI = rasListNDVI, files.names = files.names, id = 1)

id2 <- plotS2parcel(listS2parc = rasStack, listS2parcNDVI = rasListNDVI, files.names = files.names, id = 2)

id3 <- plotS2parcel(listS2parc = rasStack, listS2parcNDVI = rasListNDVI, files.names = files.names, id = 3)

---

ZADATAK

+Napraviti kratak izveštaj (RMarkodown) sa vizualizacijom i računanjem indeksa hlorofila. +Sačuvati podatke u formatu GeoTIFF za razmenu i čuvanje rasterskih prostornih podataka.